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Abstract 

Stochastic variational integrators for constrained, stochastic mechan- 
ical systems are developed in this paper. The main results of the paper 
are twofold: an equivalence is established between a stochastic Hamilton- 
Pontryagin (HP) principle in generalized coordinates and constrained co- 
ordinates via Lagrange multipliers, and variational partitioned Runge- 
Kutta (VPRK) integrators are extended to this class of systems. Among 
these integrators are first and second-order strongly convergent RATTLE- 
type integrators. We prove strong order of accuracy of the methods pro- 
vided. The paper also reviews the deterministic treatment of VPRK in- 
tegrators from the HP viewpoint. 

1 Introduction 

Since the foundational work of Bismut [1981], the field of stochastic geometric 
mechanics is emerging in response to the demand for tools to analyze the struc- 
ture of continuous and discrete mechanical systems with uncertainty [2^; [^; [l|; 



ITt |18|; [1^; [11|; [15| . Within this context the goal of this paper is to develop effi- 
cient, structure-preserving integrators for long-time simulations of constrained, 
mechanical systems perturbed by white-noise forces and torques. Our strategy 
is to employ stochastic variational integrators (SVIs) JJ. 

Variational integration theory derives integrators for mechanical systems 
from discrete variational principles (H ; 14 ; 2^; EB] ■ The theory includes discrete 



analogs of the Lagrangian, Noether's theorem, the Euler-Lagrange equations, 
and the Legendre transform. Variational integrators can readily incorporate 
holonomic constraints (e.g., via Lagrange multipliers) and non-conservative ef- 
fects (via their virtual work) p5; 16]. Altogether, this description of mechanics 
stands as a self-contained theory of mechanics akin to Hamiltonian, Lagrangian 
or Newtonian mechanics. 
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Variational integrators are symplectic, i.e., the discrete flow map they define 
exactly preserves the continuous symplectic 2-form. Using backward error anal- 
ysis one can show that symplectic integrators applied to Hamiltonian systems 
nearly preserve the energy of the continuous mechanical system for exponentially 
long periods of time and that the modified equations are also Hamiltonian [7|. 
Variational integrators are also distinguished by their ability to compute statis- 
tical properties of mechanical systems, such as in computing Poincare sections, 
the instantaneous temperature of a system, etc. 

Stochastic variational integrators are an extension of variational integrators 
to so-called stochastic mechanical systems. These systems are simple mechani- 
cal systems subject to certain random perturbations, and have their origins in 
Bismut's foundational work Bismut showed that the flow of these stochastic 
mechanical systems extremize an action integral whose domain is the space of 
semimartingales on configuration space. Bismut's work was further enriched and 
generalized to manifolds by recent work l3; U | . Lazaro-Cami and Ortega show 
that this general class of stochastic Hamiltonian systems on manifolds extrem- 
izes a stochastic action defined on the space of manifold-valued semimartingales 
Moreover, it has been shown that for a subclass of these systems, one can 
prove a converse, namely, a.s. a curve satisfies so-called stochastic Hamilton's 
equations if and only if it extremizes a stochastic action 

With this variational principle, one can design SVIs [3| • Like their determin- 
istic counterparts, these methods have the advantage that they are symplectic, 
and in the presence of symmetry, satisfy a discrete Noether's theorem. More- 
over, symplectic methods for stochastic mechanical systems on vector spaces 
have been shown to capture the correct energy behavior even in the presence 
of dissipation [13; [3]. In particular, the energy injected or dissipated by the 
symplectic integrator is not an artifical function of the timestep. Moreover, the 
energy behavior is accurately captured by the integrator. 

These structure-preserving properties are quite crucial. Consider, for in- 
stance, simulation of a simple mechanical system subject to random forces with 
amplitude e. Suppose further the unperturbed system preserves energy, mo- 
mentum, and possesses a first integral. Consider simulating this system with 
a higher-order accurate method, a standard integrator with simultaneous pro- 
jection onto energy, momentum, and first integral level sets, and a stochastic 
variational integrator. If e ^ 0, a stochastic (physical) perturbation in the en- 
ergy, momentum, and first integrals will appear. However, it is not clear how to 
modify the projection-based method to correctly capture these stochastic per- 
turbations when e 7^ 0. Moreover, the higher-order accurate method requires 
a time-step smaller than the amplitude of the perturbation in order to accu- 
rately represent its effects. And even then, if the time-span of integration is 
long enough systematic, artificial drift in these quantities will appear. 

On the other hand, the stochastic perturbation in energy and momentum 
captured by the stochastic variational integrator is mechanical, i.e., it is only 
due to the e-random forces. This is because these schemes define flows of dis- 
crete mechanical systems. Moreover, even when the amplitude of perturbation 
is not small, due to symplccticity we conjecture that SVIs on manifolds not 
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only possess good long-time energy behavior, but also perform well in comput- 
ing statistical properties, such as autocorrelation functions and the empirical 
distribution. We will investigate these questions in future work. 

As far as we can tell, the extension of these structure-preserving integra- 
tors to stochastic mechanical systems with holonomic constraints has not been 
completed. The main goal of this paper is to extend SVIs to holonomic con- 
straints, and in particular, introduce constrained, stochastic variational 
partitioned Runge-Kutta (VPRK) methods for such systems. Within this 
family we exhibit a first-order strongly convergent, symplectic integrator for 
constrained mechanical systems. We use a technique due to Vanden-Eijnden 
and Cicotti [2006] to prove order of accuracy of the integrators that appear in 



Future work will consider extensions of these schemes to Langevin equations 
with holonomic constraints. Continuous Langevin processes have been general- 
ized to submanifolds S C and shown to be ergodic The generalization 
to constraint submanifolds was done by appending holonomic constraints to 
Langevin equations. One can then check that the infinitesimal generator of the 
constrained Langevin process commutes with the Gibbsian density restricted to 
E to determine that the restricted Gibbsian measure is an invariant measure 
of the constrained Langevin process. To prove this measure is unique one uses 
standard arguments based on nondegeneracy of the diffusion and drift vector 
fields on the momentums. On this note it would be interesting to ascertain 
ergodicity of SVIs for constrained Langevin systems. 

2 Constrained VPRK Integrators 

This section is provided to fix notation and clarify some aspects of deterministic 
constrained VPRK integrators which will be relevant in the stochastic context. 
We also provide a novel proof showing that constrained VPRK integrators can 
be derived directly from a discrete variational principle without explicitly in- 
troducing a discrete Legendre transform. For more details on unconstrained 
VPRK integrators we refer the reader to . 

The setting of this section is a real, n-dimensional vector space Q and a 
mechanical system with smooth, holonomic constraint function, g : Q ^ MJ^, 
k < n, that has a regular value at G K*^. The mechanical system's configu- 
ration space is given by the constraint submanifold: S = g~^{0) C Q. We will 
introduce Lagrange multipliers to prove an equivalence between a Hamilton- 
Pontryagin (HP) variational principle on S and a constrained HP variational 
principle on Q. We will then show how to discretize this system to obtain 
variational RATTLE integrators for constrained mechanical systems. 

The variational and symplectic character of VPRK integrators is discussed 
in [U [H; I3]- In what follows we will explicitly use the HP perspective, and 
specifically, extend the results in Hairer et al. [2006] to constrained systems using 
the HP perspective. It is worth mentioning, the Hamiltonian counterparts of 
the constrained VPRK methods, so-called symplectic partitioned Runge-Kutta 
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methods, are also well understood 
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Discretization of HP Action We will adopt a HP viewpoint to describe the 
action integral of this mechanical system with constraints. The HP description 
unifies the Hamiltonian and Lagrangian descriptions of a mechanical system 
[26; 27; 2; 3]. 

Definition 2.1. The Pontryagin bundle of a manifold M is defined as PM = 
TM ®T* M . Fixing the interval [a, 6] C M and xi,X2 G M , define the HP path 
space on PM as 

CiPM,XuX2) ^ {{q,v,p) e C^{[a.b],PM) \ qia) ^ x,, = X2}. 



The HP path space is a smooth infinite-dimensional manifold. One can 
show that its tangent space at c = {q,v,p) € C{PQ, xi,X2) consist of maps 
w= {q,v,p,Sq,dv,dp) G C°°([a, 5], r(PQ)) such that 6q{a) = 6q{b) = 0. 

Definition 2.2. Fix qi,q2 G N. Define the unconstrained HP action © : 

C(PQ,gi,g2) as 



L{q,v)dt+ {P,^ 



V ) dt 



To discretize (5 we first discretize the kinematic constraint: dq/dt — v. An 
s-stage RK method is employed to discretize the kinematic constraint. Let [a, b] 
and N be given and define the fixed step size h = {b — a) /N and ti^ = hk, 
k — 0, N . The reason for using an s-stage RK discretization of the kinematic 
constraint is that the theory on such methods (order conditions, stability, and 
implementation) is mature. See, for instance, [8|. 

Definition 2.3. Consider the first order differential equation 

q^f{t,q), qiO)^qo, q{t)eQ. (1) 

Let bi, Oij G K (i, j — 1, ■ ■ ■ ,s) and let Ci — ^ij- •s-s^affe RK approx- 

imation is given by 

Ql = qk + hY,j^^^a,jf{tk + Cjh,Ql), i = l,---,s, 
Qk+i = qk + hJ2j=ibjf{tk + Cjh,Ql). 

The vectors q^ and Ql. are called external and internal stage vectors, respectively. 

It follows that an s-stage RK method is fully determined by its a-matrix and 
6-vector which are typically displayed using the so-called Butcher tableau: 
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Suppose that v{t), t e [a, b], is given. Then an s-stage RK approximant apphed 
to q = v{t) yields: 



Qk = Qk + hJ2''j=i<^'ij'>^{tk + Cjh), i = l,---,s, 
Qk+i = qk + hYfj^ibjv{tk+Cjh), A; = 0, • • • , A^. 



(3) 



In what foUows v(tk + Oih) for i — 1, s will be introduced as internal stage 
unknowns and will be determined as a critical point of a discrete action sum. 



VPRK Integrator for Constrained Systems The VPRK method will be 
derived from a discretization of the HP action integral in which the kinematic 
constraint over the kth-time step is replaced with its discrete approximant: ([3]), 
and the integral of the Lagrangian over the kth-time step is approximated by 
the following quadrature: 



tk + h 



Liq,v)dt^hbMQlV^)- 



The constraint g{q) = is enforced for all internal stage positions {Q].} using 
Lagrange multipliers as follows. 

Definition 2.4. Fix two points qi andq2 on Q and define the discrete VPRK 
path space as: 

(7(0) = qi,q{tN) = 92}, 



id the discrete constrained VPRK action sum 0^/ : 



R by: 



N-l s 



k=0 1=1 



bMQl Vk) + {pI (Ql - qk)/h - J2 <^^jVk 
\ j=i 



+ (pk+i, {qk+l - qk)/h -YJ'.VI^ + h {Al,g{Ql)) 



The following condition on the coefficients of the s-stage RK method will be 
important in obtaining a well-defined discrete update on phase space that also 
respects the constraints [lo| . 
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Condition 2.1. Consider an s-stage RK method with given b-vector and a- 
matrix. Assume bk ^ and set Okj = bj — ajkbj/bk- The coefficients of the 
s-stage RK method satisfy: 

aij = 0, flsi = 6i, foj / 0, i = l,...,s 
j ttikCikj ] invertible. 

\k=l ) ij=2 

Under this condition on the coefRcients of the s-stage RK method, one can 
prove the following theorem. 

Theorem 2.5. Given an s-stage RK method with b-vector and a-matrix that 
satisfy condition \2.1[ a Lagrangian system with smooth Lagrangian L : TQ — > M 
such that d^L/dv^ is invertible, and smooth holonomic constraint function g : 
Q — > R*^. A discrete curve Cd G Cd{qi, (72) satisfies the following VPRK method: 

Qk = qk + hJ2j=i''v'^i' 
Qk+i = qk + hY^'.^^bjVj!, 

PI - Pk + hEU{b,-'-^)[^iQlv^) + 'f^iQir.Ai), 
. 9iQl) - 0. 

(4) 

for i = 1, • • • , s and k — I, ■ ■ ■ , N — 1, if and only if it is a critical point of 
the function Q5d '■ Cd ^ R, that is, d<&d{cd) = 0. Moreover, there exist {A\}i=i 
such that the discrete flow map defined by the above scheme, : T*S T*S , 
preserves the canonical symplectic form on T*S . 

Proof. Under the assumptions of the theorem, the existence of a numerical 
solution of ^ and a discrete flow map Fh : T*S — > T*S is guaranteed. That is, 
given {qk,Pk) £ T*S one can solve dH) for {qk+i,Pk+i) <= T*S. In particular, the 
condition that d^L/dv"^ is invertible, implies one can eliminate {V^Yi=i using 
the Legendre transform of L. The condition 12.11 implies that one can determine 
{A^Li so that {qk+uPk+i) e T*5 0. 

The differential of ©d(cd) in the direction z = {{5qk,5pk},{SQ\,5Vl,piS\l=i,{5h\}l^i) 
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is given by: 



N-l s 



Pfe, Wk - 5qk)/h -^a^jSV^ ) + (pk+l, {Sqk+l - Sqk)/h -"^bjSV^ 



i=i 



hbi 



SpliQi-qk)/h-J2a,,V^j + (^ 
{6Alg{Ql)) + (^Al^{Ql).5Ql 



Spk+i, {qk+i - qk)/h -'^bjV^ 



Collecting terms with the same variations and summation by parts using the 
boundary conditions 6qQ = Sqpj = gives, 



N-l s 



fc=l i=l 



dq 



dL 



-Pk 



Pk+i +Pk-^Pk] ■ hk + h I b, — {Ql,Vk) - aj,p{ - b,pk+i ) • SV,^ 



h ( Spl, {Ql - qk)/h ~ a.jVl )+h( Spk+i, {qk+i - qk)/h - ^ bjV^ 



hb,{6A^^,g{Ql)) 



Since d&d{cd) = if and only if d©^ • z = for all z e T^^Cd, one arrives at the 
desired equations with the elimination of p\. and the introduction of the internal 
stage variables PI = dL/dv{Q\, V^) for i = 1, • • • , s. Conversely, if Cd satisfies 
Q then d<S>d{cd) = 0. 

We will employ the variational proof of symplecticity to prove that Fh is 
symplectic. This proof is standard, however, we provide it here to emphasize 
that the symplectic form that is exactly preserved by the method is the canonical 
one on T*S. 

Consider the subset of Cd given by solutions of ([4]). Let 6^ denote the 
restriction of &d to this space. Since each of these solutions is determined 
by an initial point on T*S', one can identify this space with T*S, and hence, 
&d '■ T*S R. Since ^d is restricted to solution space, 

d0d(go,Po) ■ {5qQ,5po) = {pN,SqN) - {po,Sqo) . 

Observe that these boundary terms are the canonical one-forms on r*^ eval- 
uated at fc = and k — N . Preservation of the symplectic form follows from 
d^Cd = 0. □ 
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Variational RATTLE for Constrained Systems The variational RAT- 
TLE integrator is tlie Lagrangian analog of the RATTLE algorithm originally 
proposed as a constrained version of Verlet in [2lj . It was shown to be sym- 
plectic in 13[ , and was extended to general constrained Hamiltonian systems by 
[10| . It is defined by the following two-stage RK discretization of the kinematic 
constraint (implicit trapezoidal rule) , 




1/2 




1/2 



1/2 1/2 

Given h and {(ikiPk)i the method determines {qk+iTPk+i) and two Lagrange 
multipliers Aj., A^ by solving the following system of equations, 



qk+i 

^k 



Pk+1 




Pk+1 



9{qk 



k+l, 



-I) 

ff(gfe+i) -Vk+i, 
^{qk+i,Vk+i), 



(5) 



By theorem 12.51 variational RATTLE defines a symplectic scheme. Moreover, 
as is well-known it is second-order accurate. 



Variational Euler for Constrained Systems One can relax the condi- 
tions assumed in theorem 12.51 on the coefficients of the s-stage RK method Q . 
For instance, if 5s = in the s-stage RK method one can still obtain a well- 
defined variational integrator. Consider as an example the following two-stage 
RK method, 


_ J 0_ 

1 

The corresponding variational integrator is given by: 

qk+i = qk + hvk+i, 

Pk+1 = Pk + h(^^{qk,Vk+l) + ^iqkri^l) 

= g{qk+i), 

^ Pk+1 = ^{qk,Vk+l)- 

However, the corresponding discrete flow does not satisfy the "hidden" velocity 
constraint, and hence, does not define a map from T*S to T*S. To satisfy the 
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hidden constraint a projection step {qk^i,pk+i) i-^ {qk+i,Pk+i) is taken: 



r Pk+1 = Pk+i + h^^iqk+iTAl 

I = ^{q,+,)-Vk+i, (7) 
[ Pk+1 ^{qk+uVk+i). 

One can check that this projection step defines a symplectic map. Thus, the 
composite map {qk,Pk) ^ (9fc+i,Pfe+i) ^ iqk+i,Pk+i) is symplectic. It is the 
Lagrangian version of the constrained symplectic Euler method 7]. 

Another example relaxing the assumptions in theorem 12.51 is given by re- 
garding the 2-stage RK method as being single-stage implicit Euler, 





1 




1 



The corresponding variational integrator is: 

Qk+i = qk + hvk+i, 

Pk+1 - Pk + h[^{qk,Vk) + ^^iqkr^l) 

= g{qk+i), 
, Pk+1 = ^iqk+i,vk+i)- 

We conclude with a theorem summarizing the structure-preserving properties 
of these constrained variational Euler methods. 

Theorem 2.6. The composition of one step of or (0] and the projection 
step defines a symplectic integrator on T*S. Moreover, these integrators 
are first- order accurate. 



3 Constrained, Stochastic Mechanical Systems 



Constrained, Stochastic HP Principle The setting in this section is an 
n-manifold Q and a stochastic mechanical system with smooth, holonomic con- 
straint function, g : Q BJ' , k < n, that has a regular value G M*^. Its 
configuration space is given by the constraint submanifold: S = g~^{0). In this 
section we will introduce Lagrange multipliers to prove an equivalence between 
a stochastic variational principle on S and a constrained, stochastic variational 
principle on Q. 

We will adopt a HP viewpoint to describe this mechanical system with ran- 
dom perturbations and will refer to this perturbed system as a stochastic me- 
chanical system. The HP principle unifies the Hamiltonian and Lagrangian 
descriptions of a mechanical system [2^ [23 S [1]. Roughly speaking, in the 
stochastic context it states the following critical point condition on PS (cf. def- 
inition [2lll), 



L{q, v)dt + ^ 7,(g) odW, + { p, 

i=l 



dq 
'di 



dt 
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where {q{t), v{t),p{t)) £ PS are varied arbitrarily and independently with end- 
point conditions q{a) and q{b) fixed, and : Q —f M. and Wi for z = 1, m are 
stochastic potentials and Wiener processes respectively. This principle builds 
in a Legendre transform, stochastic Hamilton's equations and stochastic Euler- 
Lagrange equations. The action integral in the above principle, consists of two 
Lebesgue integrals with respect to the Lebesgue measure dt and m Stratonovich 
stochastic integrals. This action is random, i.e., for every sample point uj € fl 
one will obtain a different, time-dependent Lagrangian system. However, each 
system possesses a variational structure as made precise in (3]. The following 
definitions will be useful to state the constrained, variational principle of HP 
for mechanical systems with holonomic constraints. 

Definition 3.1. Let i : S ^ Q be the inclusion mapping. Fixing the interval 
[a,b], define constrained HP path space as 



Cf^([a,6],gi,92) 



--{iq,v,p) : [a,b]^PS \ 
qeCH[a,b],S), iv,p)eC'{[a,b]), q{a) 



qii q{b) = q2}, 



and the unconstrained HP path space as 



C^^([a, 6], 91,92) ={{q,v,p) : [a,b] ^ PQ \ 

qeC\[a,b],Q), {v,p)eC'^i[a,b]), q{a)=^iiqi), zfe)}- 

Let P) be a probability space, J-t, t ^ [a, b], be a nondecreasing family 

of (T-subalgebras of J-", w G fl, and {Wi{t,uj),J-t), « = 1, ...,m, be independent, 
real-valued Wiener processes. In terms of these Wiener processes, we define the 
following. 

Definition 3.2. Set = L\j,g and 7^ = 7j|^ for j — l,...,m. Moreover 
define the unconstrained action 25 : fl x C([a, 6], 91, 172) ^ K as 



© = 



L{q,v)dt + J2ljil)°dWj 



dq 



dt 



and the constrained action as i3c — © 



CHP- 



The unconstrained HP path space is a smooth infinite-dimensional manifold. 
One can show that its tangent space at c = {q, v,p) G C^^{[a, b],qi, 92) consist 
of maps w — {q, v,p, Sq, Sv, Sp) E C"([a, b],T{PQ)) such that Sq{a) — Sq{b) = 
and q, Sq are of class . 

As opposed to using generalized coordinates on PS, we wish to describe 
the mechanical system using constrained coordinates on PQ and introduce La- 
grange multipliers to enforce the constraint. However, because of the stochastic 
component of the action, the standard Lagrange multiplier theorem will not 
apply directly and one cannot introduce Lagrange multipliers in the standard 
way. Instead, we will introduce the Lagrange multiplier using the following. 
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Definition 3.3. Given /i e C°{[a,b],T*Q) and /2 e C\[a,b],TQ), define 

Differentiability of the flow map on T* S will be defined in the mean-squared 

sense. In the following wc define mean-squared derivatives on a Banach space E 
with the understanding that this notion can be extended to any manifold using 
a local representative of the flow map. 

Definition 3.4 (Mean-squared Derivative). The mean squared norm of f : 
E X CI E is given by: 

\\f{x,u;)\\ = {E{\f{x,u;)f)f' 

Using this norm one can define the derivative of f in the standard way, i.e., 
f is m.ean squared dijferentiable at a G E if there is a bounded linear map 
Df(a,uj) : E E that satisfies, 

,. \\f{a + 6,w)-f{a,u;)-Df{a,u;)-6\\ 



With the above definitions we prove the following. 

Theorem 3.5 (Constrained, Stochastic HP Principle). Given a stochastic 
mechanical system with Lagrangian L : TQ — » M such that d'^L/dv^ is invertible, 
stochastic potentials 7j : Q ^ R for j = l,...,m, and holonomic constraint 
g : Q ^ M.'' with S = g~^{0). Then the following are equivalent: 

(i) z = {q,v,p) e C^^{[a,h],qi,q2) extremizes <&c- 

(ii) z = {q,v,p) € C^^{[a,b],qi,q2) satisfies stochastic HP equations 
dq = vdt, 

dp = 0L^(q,v)dt + j:7=l'ilil)°dWj, (9) 
dv 



(iii) There exists \ G C°([a, 6], K'') such that z = {q,v,p) E C^^ {{a,b],i{qi),i{q2)) 
and A extremize the augmented action ©(z,A) = ®(z) + (A, $(2)) where 
^{z){t) = dg/dt{q{t)) and (A, <i>{z)) = j'^ {X{s), $(z)(s)) ds. 

(iv) There exists Xe C°([a,6],M'=) such that z = {q,v,p) G C"^{[a,b],i{qi),i{q2)) 
and A satisfy the constrained, stochastic HP equations 



dq = vdt, 

OL, . (10) 



dp = ^iq,v)dt + Y:U^{l)°dWj + ^^{qr -dX, 

p = ^M, 



[ 9{q) = 0. 
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Moreover, the flows of 0) and UU\) are mean-square symplectic. 



Remark 3.1. The constrained, stochastic HP equations should be thought of 
in integral form. In particular, the Lagrange multipler term in ilU\) by defini- 
tion \3.3\ satisfies: 



la \9q 
for any w G M'' . 



(g)* ■dX,w) ^ (A, ^{q)w 



Proof. The stochastic HP principle states that (i) and (ii) are equivalent 
Assume that (iii) is true. Then for all Vz = {vq,Vv,Vp) £ TzC^^, 

d(S(z, A) • ivz,vx) = d6(z) • Vz + {vx,Hz)) + (A, d$(z) ■ Vz) = (11) 

Since $(z) = for z e ^ and Tq^t)9{q{t)) ■ Vq{t) = for Vq{t) e TqS, 

d0(z, A) • ivz,vx) = d0{z) -Vz^O, y vze TzC^^ 

which implies (i) and (ii). 

Moreover, expanding ([TT|) and setting {vq,Vy,Vp) — {Sq,Sv,Sp) and v\ — 5X 
yields, 



dC5(z, A) • {5q,6v,6p,SX) = 



dp, ^ - i!^ ds + (^p, - Sv'^ ds + (^SX, ^gil)^ ds 



.(.|(|..,l).. 



Consider the terms involving Sp, Sv and SX. Since these variations are arbitrary, 
the following holdj^ 



dq dL d 



Since g{q{a)) — g{i{qi)) — 0, d/dt{g{q)) — imphes that g{q{t)) — for all 
t e [a,b]. 

Collecting the variations with respect to Sq in the differential gives. 



dL 



,dq 



dt 



3 = i 



dq 



d (dg 



dt \dq 



= 0. 



Using definition 13.31 we introduce the following function / : C^^{qi, 92, [a, b]) x 
C{[a,b],R'') X C^{[a,b],TQ) ^R, 



i{q,v,p,x,f) 



'-^ds + ±'-2LodW,-dp-'/ ■dX\.f 



dq 



dq 



^This follows from the basic lemma that if /, g G C''([a, b],M) and g is arbitrary then 
flf(i)g{t)dt = ^ f{t) = V tela,b]. 
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In the following it is shown that if /(g, v,p, A, /) = for arbitrary / of class 
then c = {q,v,p,X) satisfies (fTO)) . 

Let {Ua,ga} be a partition of unity on PQ x R'^. Expand / in terms of this 
partition of unity, 



dq 



Since the curves z = {q,v,p) and A are compactly supported, only a finite 
number of the ga are nonzero. For each ga nonzero, the terms in the integral can 
be expressed in local coordinates. Observe that since dq = vdt, the Stratonovish- 
Ito conversion formula implies that. 



dq 



Sq o dWj 



dq 



■ SqdWj 



for j = 1, m. 

We will select / to single out the /?th-component of the covector field in /. 
Introduce the following function /i : K ^ M for this purpose: 

t 

Observe that /i(0) = 0, h{e) — 1, and ft.'(e) = 0. Let {ep} be a basis for the 

model space of Q. Now fix /3, and define e C^([a, h]^ TQ) in local coordinates 
as: 

h{s — a)ep if a < s < a + e, 

ep if a + e < s < t — e, 

h{t — s)e^ if t — e < s < t, 

if t < s < 6. 



Ms) 



Introduce the following label to simplify subsequent calculations 

A{S): 



^{q{s),v{s))ds + J2 dW^is) - Ms) - ^ilis))* ■ dX{s)]-ep. 



In terms of one can write 

ra-\-e 



h{s-a)gais)A{s)+ I ga{s)A{s) 



i{q,v,p,x,fe) = ^ 

a 

We will show in the mean squared norm, 

lun/(g,u,p,A,/e) = ^ / gaA{s) =: I* . 



h{t - s)ga,{s)A{s) 



(12) 
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Using this result and the Borel-Cantelh lemma, one can deduce there exists {e„} 
that converges to such that I{q, v,p, A, a.s. converges to I* . It follows that 
/* = almost surely. 

We proceed to prove Since (a + 6)^ < 2a'^ + 26^, 

2 



< 2 



E 
E 



(1 - /i(s - a))5aA(s) 



(1 - h{t - .s))g^A{.s) 



(1 - h{s - a))gaA{s) 



We will only show how to bound the first term since bounding the second term 
is very similar. By continuity of {q, v,p, A), one can pick e > small enough so 
that the support of (g, v,p, A) lies in a single chart. Therefore, 

2 



E 



a+e 



a+e 




< 4 



(1 - his - a)) 



(1 — h{s — a))-^^^ds 



4E 



9(7 



a+e 



(1 — h{s — a))c?p • 6/3 



+ 4 



a+e 



(1 - h{s - a)) 



\dq 



■ dX, ep 



Since is continuous on s e [a, a + e], the first term can be bounded. 



(l_/i(s_a))_ds 



< 



Similarly, by the Ito isomctry and since is continuous on s G [a, a + e], the 
second m terms can similarly be boundeS, e.g., the jth Stratonovich integral 



can be bounded as follows, 



{1 - h{s ~ a))^^' 









dW, 




r 









{l-h{s-a)) 



dqP 



ds < 



Using (|3.3p and the integral mean value theorem, there exists a real constant 
Ci S [0, 1] such that: 

2 



(1 — h{s — a))dp ■ 6/3 



(1 - h{s - a))pp{s)C' + / pp{s)h'{s ~ a)ds 



= \\-pp{a) +pp{a + cie) 
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Similarly, there exist constants C2, C3 G [0, 1] such that 

2 



a+e I Qg* 

{\-h(s-a))(— ■d\,ei3 



A, ^ • h'{s - a)d. 



d dg 

(1 - - a))Xp — -^{q{s))ds 



dg 
dq 

dg. 



(q(a))* • X{a),ep 



£ X / \ 9g , , 

-A/5(a + C2e)^(g(a + C2e)) 



+ ^^(9(0 + cse))* • A(a + C3e),e0 

Since p/3 and A are of class C'^ and g is smooth, as e — > these terms vanish. 
Since (3 is arbitrary we have proved (|12p . Therefore, almost surely: if c = (z. A) 
is a critical point of & then c satisfies p^ . 

On the other hand, assume (iv) is true. Then almost surely: if c = (z, A) 
satisfies the constrained, stochastic HP equations, then it is a critical point of &. 
This direction is easy to confirm, since as a solution to the constrained, stochas- 
tic HP equations c is a measureable diffusion process. In fact, this direction 
is similar to the one Bismut originally established, namely that the solution of 
stochastic Hamilton's equations extremize an action function; albeit a different 
stochastic action is used in this proof [H . 

Assume that (i) is true, i.e.. 

Define = {g £ C\[a,hlQ) \ q{a) = i{qi), q{h) = ^(g2)}. Let = 
g{q{t)) where : C^([a, 5], R'"'). By the local onto theorem, there exist 

charts ([/, Lp) of and a open set F C M*' such that o Lp-^ : U' x V ^ V is 
a projection onto the second factor and (p{q) = [x, y) £ U' x V for q eU. Set 
Vq = Tx^yLp^^{vx,Vy) where € U' and Vy G V. Set ©(g) — &{q,v,p). Then, 



d©(z. A) ■iv,,vx)^T0{q)-v, + {\, T^iq) ■ -{v,) 



T2i(&op-'){x,y)-Vy + {\,-ivy) 



However, there exists A such that, 

T2{eop~^){x,y)-Vy + (^X,j^{vy) 
for all Vy V . To see this expand this expression. 



^{x,y),v),Vy'^ dt + Y^^(^^{ip ^{x,y)),v. 



V I ° dWj 



= 
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and note that a solution is given by: 

d\ = V),v)dt + ^ ^(^"'(^, V)) ° dW, 

which is an SDE for A, with globally Lipschitz drift and diffusion vector fields. 
Hence, there exists A such that, 

d©(z,A)-(i;„i;A) = 0, V e T,C^^, i;;, e rAC^°([a, 6],R^). 

Mean-square symplecticity of the flows of stochastic HP equations, and 
hence, constrained stochastic HP equations (by the equivalence) is established 
in d. □ 

Eliminating the Lagrange Multiplier One can eliminate the Lagrange 
multiplier in (fTO|) by taking the second Stratonovich differential of the constraint 
5(9) = 0: 

dg , . dg 

and, 

By differentiating the Legendre transform in (fTO)) we obtain, 

dv - (—\ ' dp - (—] ' (-^] vdt 
\d^v J \d'^v J \dqdv J 

Substituting dp from pU|) and the above into gives, 

^fyy)dt + ^ (—y^(—dt + f"^dW+^*dx\ 
d'^q ' dq \ d'^v J \ dq ^ dq ^ dq I 



dg fd^L\ ( d^L 
dq \ d'^v 



vdt = 0. 



Set P(g, v) (^^{q, v)j ^il)* and solve for dX in the above to obtain: 



dq \ d'^v J \dqdv 
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Set B{q,v) - ||(g)*P(q,^,)-i||(g) Substitute dX into m gives, 

dq = wrft, 



dp = (Id-B(g,r;))(f(g,t;)di + E;Li^(9)°'^W^., 

-f|(g)*P(g,«)-i§:f(,;,z;)dt + B(g,«)(^(g,«))t;di, 



^ ~ at, 

Remark 3.2. One can check using properties o/B and P, and the Stratonovich- 
Ito conversion formula that the Stratonovich correction term in vanishes, 
and hence, |_?^[ ) is equivalent to: 



dq — vdt, 

dp = (m-Biq,v))[^iq,v)dt + j:';^,^{q)dW, 

-^^(q)*F(q,v)-'^^iv,v)dt + Biq,v)[^Xl^v))vdt, 
dLi 



(15) 



4 Constrained, Stochastic VPRK Integrators 

General Case Let [a, h] and N be given, and define the fixed step size h = [b— 
a)/N and tk = hk, k = 0, N. The stochastic VPRK method will be derived 
by discretizing the kinematic constraint in the stochastic HP action integral 
with an s-stage RK approximant (cf. ([3])), and the integral of the Lagrangian 
over the kth-time step by the following quadrature: 

L{q,v)dt^hbMQl:V^)- 



Motivated by the methods introduced in [18|; [l9|, the stochastic part of the 
stochastic HP action will be approximated by: 

rtk+h 

-friq) O dW « {Vi(f>r + Killer) IriQl), 

'tk 

where (j>r : ^ ^ M. and 4'r ■ ^ R for r = 1, ...,m. With this discretization 
3/2-order strongly convergent methods can be derived. The constraint g{q) = 
is enforced for all internal stage positions {Q^jf^i using Lagrange multipliers 
as in the deterministic case. 

Definition 4.1. Fix two points qi and q2 on Q and define the discrete con- 
strained, stochastic VPRK action sum (&d : Q, x Cd ^ R by: 

N-l s 



fc=0 i=l 



hbMQl V^) + {Mr + ^^^4>r) 7r(gl) + kb, {A^ g{Ql)) 



-h [pI, {Qi - qk)/h - J2 ""^i^kj + ^ {Pk+1: {qk+1 - qk)/h - J2 ^J^s 
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In the following we assume that the VPRK method has a well-defined solu- 
tion much like the assumption made in theorem 12.51 



Theorem 4.2. Given an s-stage RK method with b-vector and a-matrix, a 
Lagrangian system with smooth Lagrangian L : TQ — )■ M such that d^L/dv^ 
is invertible, stochastic potential 7^ ; Q — > K for r = l,...,m, and smooth 
holonomic constraint function g '■ Q ~^ M'^. A discrete curve Cd G Cd{qi,q2) 
satisfies the following stochastic VPRK method: 



Ql 
qk+i 

pi 



Pk+1 



qk + hY.'j^i a-tjV^ , 

qu + hTU^oVl. 

vk + hj:U - ^) {^(Qlvi) + %{QiY ■ Ai) 



Vk + h [waiQiyi) + %iQiy ■ K 



ELi E^i i'^j'Pr + Kji^r) ^{Qi), 



PI = 



9{Ql) 



0. 



(16) 

for i — 1, • • • , s and k — 1, ■ ■ ■ , N — 1, if and only if it is a critical point of 
the function : Cd — > M, that is, d(3d{cd) = 0. Moreover, there exist {A^}|^j^ 
such that the discrete flow map defined by the above scheme, Fh : T* S ^ T*S, 
preserves the canonical symplectic form on T*S. 



Proof. We assume the existence of a numerical solution of (jl6l) and a discrete 
flow map Fh : T*S -> T*S. 

The differential of ©d(a;, Cd) in the direction z {{Sqk, Spk}, {SQl, SV,^,pl}f^i, {dA^jf^^] 
is given by: 



d6d(w,Cd) • z 

N-l s 

Y^Y^hb^ 



k=0 i=l 



FIT BT 

— {Ql, Vl) ■ SQl + g^iQiX) ■ SV^ 



Pk 



{6QI - Sqk)/h - X^a^j'^^fe ) + (Pfc+l' i^lk+l - Sqk)/h -Y^J^^i 



{Ql - qk)/h - Y ^^j^k ) + ( ^Pk+i, {qk+i - qk)/h ~ Y bjV^ 



hbi 



{SAl,g{Ql)) + {Al,^{Ql)-SQl 



Collecting terms with the same variations and summation by parts using the 
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boundary conditions dqo = Sqn = gives, 



de5d(w,Cd) • z 

N-l s m 



k—1 i—1 r—1 



dv 

1=1 / \ j=i 



+ h (^6pl, {Ql - qk)lh -Y/'.jVl'j + h {^Spk+i, {qk+l - qk)/h - ^ b^Vi j 
+ hb,{SAl,g{Ql)). 

Since d(8d{i^,Cd) = if and only if d&d ■ z ^ for all z G Tc^Cd, one arrives 
at the desired equations with the elimination of p\, and the introduction of the 
internal stage variables PI = dL/dv{Q\, T^-) for i = 1, • • • , s- Conversely, if Cd 
satisfies (fT6)) then d0d(w,Cd(a;)) = 0. 

The proof of symplecticity using the variational principle is very similar to 
the proof of theorem 12.51 and is therefore omitted. □ 



Constrained, Stochastic Variational Euler As an example we consider 
the following simple first-order strongly convergent integrator. One can derive 
higher-order accurate integrators for constrained systems following the proce- 
dure for unconstrained stochastic Hamiltonian systems described in 



tne proc 

US. 



Let ~ A/'(0, h) be normally distributed random variables for r — 1, 



, m 



and k = 0,...,N — 1. The constrained, stochastic variational Euler inte- 
grator is given by: 

qk+i ^qk + hvk+1, (17) 
Pk+i =Pk + h^{qk,vk) + J2 + h^^iqkTAl (18) 

0=5(9fe+i), (19) 
Pfe+i =-^('7fc+i,Wfe+i). (20) 

together with the projection step Given {qk^Pk) and ft,, this integrator deter- 
mines {qk+iTPk+i) by eliminating Vk+i using the Legendre transform, eliminat- 
ing p^+i using (HH]), and determining by satisfying the constraint g{qk+i) = 0. 
One then takes a projection step by solving ([7]) for Pk+i- The integrator defined 
by the composite map, {qk,Pk) ^ {qk+i^Pk+i) ^ {qk+i,Pk+i), has the following 
properties. 

Theorem 4.3. Constrained, stochastic variational euler is mean-squared sym- 
plectic and first- order strongly convergent integrator for mU\) . 
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Proof. We will prove this by using the technique provided in [23| . 

Consider the following first-order strongly convergent, Euler-Maruyama in- 
tegrator applied to (Hi 



<ik+i = qk + hvk, 

Pk+i = Pfe + (Id - B(gfe, ^;fe))(/.f(9fe,^;,) -1-1:7=1 ^(<Zfc)S}- 

-h^{qk)*'Piqk,Vk)~'^^ivk,Vk) + hB{qk,Vk) (^j(afc,'ffe) • Vfcj 

. Pk+l = ^{Qk+l,Vk+l)- 

(21) 

The order of accuracy of constrained stochastic variational Euler will be deter- 
mined by checking that the update determined by a single step of (fT7l) - ((20| with 
the projection step agrees with a single step of (PT|) to 0{h?). 

For this purpose the following expansions to 0{h'^) are introduced: 

qk+i =qk + hvk+i, Vk+i = Vk + hvl^j^, pk+i ^Pk + hpl^i 
Vk+i = Vk + /iWfc+i, Pk+i ^Pk + hpl+i 

Substituting these expansions into (|17p confirms that the position update agrees 
to 0{h^). However, checking that the momentum or velocity updates agree to 
0{h?) is more involved. To do this expand in a Taylor series about g^, and 
use g{qk) = and f|(gfc) • -y^. = to obtain, 



5(gfe+i) = h'^^iqk) ■ iUi + Y|J^(ft)(«/c,t'.) + 0{h^) ^ (22) 
From (fT8|) it is clear that. 



hp\^, ^h^{qk,Vk)+p^ ^{qk)B^ + h^iqkTAl (23) 

Moreover the Taylor series expansion of (PO)) gives, 

h^iqk,Vk) ■ Wfe+i + h——{qk,Vk) ■ Vk = hpl^^ + Oih^) (24) 

Using invertibility of §3^(gfe,Wfc), one can rewrite ((24|) as, 

hvl+i = -h—r-[qk,Vky^——{qk,Vk) ■ Ufc + h—^{qk,Vky^pl_^,^ + 0{h^) 
O'^v oqov O'^v 

Substitution of (|23|) into the above equation gives. 



d'^v J I, \ dqdv 



k 
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Substituting the above expression into gives, 



l)..(S)r(''f<-"''\|:?7<""''-''l'"'-*' 

From which it follows that the Lagrange multiplier satisfies, 



Substitution into (f23|) gives 



To determine p^,).]^ we expand the hidden velocity constraint in ^ to obtain, 

^{qk+l)vk+l = h^iqk) ■ vl^, + h^^(qk){vk,Vk) + O(h') = (26) 
Expanding the momentum update in ([7]) gives, 

hpU = hpl+, + h^iq.YAl + 0{e). (27) 
Expanding the Legendre transform in ([7]) gives, 

- ^ (S)r i'"'*^'-" (£)r (S). "-°""'> 

Substitution of the above into (|26p yields, 
+ /''T^(%)(«fe,i^fe) + 0(/i') = 
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Solving the above for and substituting into ([77|) gives, 
Substituting ([25]) into the above and using the following identities 



(Id-Bfc)Bfc =0, (Id-Bfc)^((jfc)*p-l^(qfe)(i;,,t;fc) = 



implies that, 

+ ^Bfc (g;)^ • - ;^|(q,)*P,-^|J^fe)K, + 0{h^) (28) 

which agrees with to ©(/i^). 

Mean-square symplecticity of the map (qk,Pk) '-^ ilk+iiPk+i) follows from 
the proof of theorem 14.21 The projection step is also mean-square symplectic. 
Therefore the composite map is mean-square symplectic with respect to the 
symplectic form on T*S since the space of symplectic maps forms a group. □ 
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Example 4.1 (Stochastically Perturbed Spherical Pendulum). Consider the 
spherical pendulum where Q = M.'^ and S — with holonomic constraint given 
by the zero level-set of g{q) = WqW^ — 1. Let {61,62,63} denote an orthonormal 
basis for R'^ with 63 corresponding to the direction of gravity. Its Lagrangian is 
given by: 

Consider the following stochastic potentials: ^i{q) = sin(q • Ci) for i = 1,2,3. 
Applying stochastic variational Euler to the resulting stochastically perturbed 
spherical pendulum one obtains the strong error plot shown in Figure\^ confirm- 
ing theorem \4.3\ Since a strong order of convergence implies at least the same 
weak order if not higher, it is not surprising that the weak error plot shown in 
Figure\^ seems to predict 3/2-order weak convergence. The dashed lines in the 
plots give appropriate reference slopes in each case. 




Figure 1: Strong error plot for stochastically perturbed spherical pendu- 
lum. A log-log graph of the strong error of stochastic variational Euler applied to 
the stochastically perturbed spherical pendulum. Green and cyan solid lines repre- 
sent the strong errors in momentum and position respectively. For reference a dashed 
line of slope 1 is included. Observe that the order of convergence is consistent with 
theorem 14.31 
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Time-Step Size 



Figure 2: Weak error plot for stochastically perturbed spherical pendu- 
lum. A log-log graph of the weak error of stochastic variational Euler applied to the 
stochastically perturbed spherical pendulum. Green and cyan solid lines represent the 
strong errors in momentum and position respectively. For reference a dashed line with 
slope 3/2 is included. This order of convergence is also consistent with theorem 14.31 
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